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ABSTRACT 

We give a full nonlinear numerical treatment of time-dependent 5d braneworld geometry, 
which is determined self-consistently by potentials for the scalar field in the bulk and at 
two orbifold branes, supplemented by boundary conditions at the branes. We describe the 
BraneCode® , an algorithm which we designed to solve the dynamical equations numerically. 
We applied the BraneCode to braneworld models and found several novel phenomena of the 
brane dynamics. 

Starting with static warped geometry with de Sitter branes, we found numerically that this 
configuration is often unstable due to a tachyonic mass of the radion during inflation. If 
the model admits other static configurations with lower values of de Sitter curvature, this 
effect causes a violent re-structuring towards them, flattening the branes, which appears as 
a lowering of the 4d effective cosmological constant. 

Braneworld dynamics can often lead to brane collisions. We found that in the presence of 
the bulk scalar field, the 5d geometry between colliding branes approaches a universal, ho- 
mogeneous, anisotropic strong gravity Kasner-like asymptotic, irrespective of the bulk/brane 
potentials. The Kasner indices of the brane directions are equal to each other but different 
from that of the extra dimension. 
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1 Introduction 



Braneworlds embedded in higher dimensions bring new powerful concepts in cosmology pQ, 
as well as they do in fundamental superstring/M theories and phenomenological high energy 
particle physics [3EHUIS1- Branes enrich our view with new ideas underlying the four di- 
mensional effective field theory, bringing, most importantly, new geometrical images beyond 
it. For instance, in this context the four dimensional cosmological constant is the curva- 
ture of the brane, and we should explain why the brane we live in is almost flat. Inflation, 
by contrast, corresponds to curved branes. There are interesting ideas for realizing early 
universe inflation in braneworld scenarios where, for example, concepts such as the inflaton 
potential and inflaton decay are reformulated in terms of brane-brane or brane-antibrane 
interactions |6 , or topics of brane collisions. 

The language and images of the braneworld theories are commonly shared by fundamental 
and phenomenological high energy physics theories, general relativity and brane cosmology, 
with different degrees of trade between dynamics and simplification. 

The compactification of the extra space is often a key issue in brane models. For example, 
a stable radion field, controlling the volume of the extra space, is usually needed to recover 
standard four dimensional cosmology at "late" times, and to fulfill precision tests of general 
relativity. In addition, the compactification has to be consistent with the fact that bulk fields 
have not yet been excited in accelerator experiments, because they are too massive and/or 
too weakly coupled to the visible brane. Schemes for compact inner dimensions typically 
rely on the interplay between bulk and brane dynamics. 

So far, the control of dynamical, time-dependent, cosmologically relevant solutions in 
the fundamental, comprehensive theory has been rather limited. Relatively simple, yet 
meaningful, are the five dimensional phenomenological braneworld models with two orbifold 
branes at the edges, where our 3 + 1 dimensional spacetime is one of the branes embedded 
in the (warped) five dimensional space. These models often include one or more bulk scalar 
field(s) <p with the potential V(4>) and self-interaction potentials Ui(<p) at the two branes, 
as well as other fields \ localized at the branes. This class of braneworld models covers 
many interesting constructions including the Hofava-Witten theory the Randall-Sundrum 
model with a phenomenological stabilization of branes [7||H], warped geometry with bulk 
scalars [1011201, supergravity with domain walls j^], and others. 

There is a number of important papers studying static geometries with branes, including 
flat stabilized branes, in agreement with low energy physics, curved de Sitter branes, corre- 
sponding to early universe inflation, and small fluctuations around static warped geometries. 
The cosmological evolution has been studied in some of these pioneering works in the sim- 
plest cases in the absence of any scalar field [T21ll3| . The 4d evolution on the brane, in terms 
of effective Friedmann equations, is typically different from the standard four dimensional 
cosmology. The effective 4d Einstein equations on the brane were also derived for the more 
general situation of self-consistent geometry with the bulk/brane scalar field [Tlj . 
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Standard cosmology can be recovered after the extra dimensions have been stabilized pTj. 
In this respect, the presence of bulk scalar field(s) becomes crucial. In this more relevant 
case, however, the evolution is only known for limiting situations. In general, the system is 
very complicated, since the effective four dimensional Einstein equations are not closed and 
require solutions of the full five dimensional equations [31] . 

In this paper we address the problem of self-consistent fully nonlinear dynamics of the 5d 
braneworld with bulk scalar field(s) with a bulk potential as well as brane potentials required 
for brane stabilization. We consider plane-parallel orbifold branes, so that the problem is 
effectively two dimensional, with the metric and the fields depending on time and the extra 
dimension (t, y). Although this setting is already too involved to be studied analytically, 
it is still significantly simpler than what we will need in order to understand cosmological 
solutions in "realistic" higher dimensional theories. However, as we shall see, already this 
step requires the introduction of new techniques. We have designed and used a numerical 
code to solve the partial differential equations describing the system of non-linear gravity and 
a scalar field, complementing the existing approaches to this problem found in the literature. 

We aim for generic features of braneworld dynamics, in particular, attractor solutions. 
They will generally depend on the specific braneworld model, i.e. on the bulk/brane scalar 
field potentials V(<j)),Ui((j)). As a simple illustration, consider a static five- dimensional 
warped geometry with a bulk scalar and four dimensional slices of constant curvature. It 
is possible to exhaust the global properties of static warped geometry using the method 
of phase trajectories [TJj, although some details of the phase portraits depend on the bulk 
potential. For this problem the phase space is three dimensional, the critical points (like 
attractors, repulsors, and others) can be identified, and all trajectories (solutions) start and 
end at critical points. 

The (t, y) problem of the time-dependent braneworld dynamics is much more complicated 
than the static (y) problem. Using the BraneCode we were trying to give examples of 
interesting dynamical features. We notice several novel phenomena including a transition 
between different warped states and a generic strong gravity solution of colliding branes. 

The plan of the paper is the following. 

In Section we give the setup of the braneworld models and write down the bulk 
equations supplemented by the junction conditions at the branes. We pay especially close 
attention to the choice of the gauge in order to have a suitable metric for the numerical 
calculations. It turns out that, without any loss of generality, it is possible to choose coor- 
dinates where the two branes have fixed positions along the fifth direction y. The geometry 
is described by two metric components A(t,y) and B(t,y). 

In Section |HJ we describe the BraneCode, an algorithm we use to solve the dynamical 
equations numerically. At the moment we have slightly different implementations of the 
BraneCode (in C + + and Fortran— 90) in order to cross-check them. We plan to release the 
BraneCode using the most optimized and documented version (in C + +). As is typical for 
numerical GR problems, we have to take initial conditions for the metric and fields which 
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satisfy the constraint equations at an initial time hypersurface. In Subsection 13. 3j we discuss 
how to fix the initial conditions for the numerical integration with the BraneCode. 

In Sections EHEl we apply our BraneCode to three braneworld models where we encounter 
qualitatively different dynamics. 

In order to check our numerical code, in Section 0] we first apply it to a simpler brane 
model without a scalar field, for which analytic solutions are known. As a playground here 
we use the Randall- Sundrum model of two branes embedded in an AdS 5d background. In 
Subsection 14.11 we first use the static RS solution without moving branes. In Subsection 14.21 
we extend the calculations to the case of moving branes. In this case the 5d geometry is 
described by the analytic AdS-Schwarzschild solution (with the mass of the virtual 5d black- 
hole screened by the branes). We compare our numerical calculations with the analytic 
solution. 

In Section El we consider de Sitter (inflating) branes which are initially in a static config- 
uration. It turns out that we often observe an instability of the inflating branes. Analytic 
calculations of small scalar perturbations around this background geometry show that the 
radion mass square m 2 for this case can be negative jTHj. A strong tachyonic instability 
predicted analytically is in full agreement with the instability of inflating branes found nu- 
merically. For certain configurations of potentials V(<f)), Ui(<f)), we find the existence of two 
warped geometry solutions with different values of the 4d cosmological constant A4 (i.e. 
the curvatures of the 4d slices). The brane configuration with the higher 4d curvature is 
in general unstable due to this tachyonic radion mode, and violently re-configures to the 
second static configuration with lower 4d curvature. We illustrate this effect with numerical 
simulations as well as analytic calculations, see also [T6j . 

In Section |B] we give an example where the instability of the brane configuration causes 
a brane collision. In Subsection 16.11 we show that the space-time metric of the 5d geometry 
between colliding branes becomes homogeneous, i.e. y-independent. The time- dependent 
solutions asymptotically cease to feel the scalar field potentials V(<f>) and Ui(<f>), and approach 
a universal asymptotic. It sounds natural a posteriori that this universal asymptotic is 
nothing but a Kasner-like asymptotic with a scalar field, which we describe in Subsection l6.2l 
The effect of the branes here is manifested by the fact that the Kasner indices associated with 
the three brane directions are equal, but different from that associated with the y direction. 
This is a strong gravity regime, so it is not surprising that the 4d induced metric on the brane 
is different from that derived with moduli approximations in terms of 4d effective theory. 

In the Conclusion we summarize the most interesting physical results. Technical details 
are collected in several Appendices. 



4 



2 Setup 



The class of braneworld models we are interested in is characterized by the following action 



where k\ = is a 5d gravitational constant. In this convention is measured in units of 

/% 1 and physical potentials are multiplied by k§ 2 . The first term describes gravity in the 
five dimensional bulk space. We use the "mostly positive" metric signature. The second 
term corresponds to a (minimally coupled) bulk scalar field with the potential V((f>). The 
last term corresponds to two 3 + 1 dimensional branes, which constitute the boundary of 
the five dimensional space. We allow for a potential term U(<p) for the scalar field at each 
of the two branes. We denote by 7 the induced metric on the two branes, and by K their 
extrinsic curvature. Here and in the following, [Q] = Q (y + ) — Q (y_) denotes the jump of 
any quantity Q across a brane (± defined respect to the normal of the brane). We assume 
5' 1 /Z 2 symmetry across each brane. 

The algorithm we have written is implemented for generic bulk and brane potentials. In 
this paper we specify the potentials introduced for the brane stabilization [Jj. We choose 



where fa is the value of <fi on the i— th brane. A 5d cosmological constant A in the bulk and 
tensions A, on the branes are included in the potentials. 

The two branes are assumed to be parallel. We denote by y the coordinate transverse to 
them, and by x the three spatial longitudinal coordinates. We assume isometry along three 
dimensional x slices including the branes. We have to specify a metric qab that respects 
this symmetry. In brane cosmology it is customary to use the metric in the form ds 2 = 
—n 2 dt 2 + a 2 dx. 2 + bdy 2 , where the metric components n, a, b depend on (£, y). However, this 
form of metric does not exhaust the freedom of the coordinate choices. Most significantly, in 
this metric the branes do not stay at the fixed positions; in general y^ = yi{t). There are other 
gauge choices, which were used for specific braneworld problems, for example coordinates 
comoving with one of the branes, the choice of the bulk scalar <fi as the y-hyper surf ace 
and others. In these contexts, a gauge in which the position of one of the two branes is 
time dependent was often preferred, and identified with the radion field 1Z (t) associated 
with the extra dimension. Although this choice may lead to an easier interpretation of the 
interbrane distance, the resulting bulk and junction conditions (see below) are significantly 




8=1,2 



(1) 




(2) 
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more complicated. In addition, in terms of the four dynamical quantities a, b, n, 0, and 7Z 
the system is actually under-determined, and some gauge fixing is needed to have a closed 
set of equations. 

For numerical simulations, it is preferable to have coordinates where neither brane is 
moving, although it is not obvious a priori that such a gauge can be constructed without 
loss of generality. It is possible to choose coordinates such that the bulk metric has the "2d 
conformal gauge" , 

ds 2 = e 2B ^(-dt 2 + dy 2 )+e 2A ^dx 2 , (3) 

This gauge still has the residual freedom to change (t, y) — > (t', y') in a way that preserves the 
2d conformal form. It can be demonstrated that this freedom can be used to fix the position 
of the two branes along y . Without loss of generality, we can locate them at y = 0, 1 . We 
found that in the 2d conformal gauge the set of bulk equations (0) acquires a relatively simple 
form, which is well suited for the numerical scheme we have adopted (see the next section). 
The possibility of choosing a gauge, in which the metric is 2d conformal and in which the 
branes are at a fixed position is shown explicitly in Appendix IA.11 As it is discussed in 
Appendix IA.21 even these requirements do not fix the gauge choice completely. 

Although in the system of coordinates we have chosen the branes are always at a fixed 
position along the y axis, their physical distance is encoded in the metric component B , 
which is a time dependent quantity. Clearly, the distance between two extended objects 
is not an invariant quantity in general relativity, and different definitions can be adopted 
when they are in relative motion. A simple heuristic possibility, which we adopt here, is to 
integrate the line element across the extra dimension at a fixed time 

D(t)= [ dyV5to= [ dye 3 ^ . (4) 
Jo Jo 

One can check that D (t) is invariant under the residual gauge freedom in our coordinates 
(jHJl, which is discussed in Appendix IA. 21 (but not under general coordinate transformations). 

For the output of our numerical calculations, we rely on gauge invariant quantities such 
as the invariants of the 5d Weyl tensor CabcdC abcd , the curvature scalar R and others. 
These invariants are calculated using the metric in the form (j3J). Additionally, we can use 
the 4 + 1 split of the 5d curvature, symbolically written as R = _R 4 + K 2 , where i? 4 is 
the curvature of the 4d slices. This will be especially useful when we work with de Sitter 
(inflating) branes of constant curvature. 

In the gauge we have chosen, the nontrivial five- dimensional Einstein equations can be 
split into three dynamical equations 

A - A" + 3A 2 - 3A' 2 = - e 2B V , 

3 

B-B" - 3A 2 + 3A' 2 = -¥- + f-- -e 2B V , 

2 2 3 

<j> - 0" + 3i0 - 3A'<f>' + e 2i % = , (5) 
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plus two constraint equations 

-A'A + B'A + A'B - A' = \<P4>' , 

2A' 2 - A'B' + A" - A 2 - AB = — — - — - -e 2B V . (6) 

6 6 3 v ; 

Dots and primes denote derivatives with respect to t and y , respectively. It is easy to show 
that the constraint equations are preserved by the dynamical equations. 

In addition, from the boundary terms in the action for the two branes we recover the 
following junction (Israel) conditions 

[A'] = T\ue B , 
[B'] = T~ U e B , 

W] = ±e B U^, (7) 

where the upper and lower signs refer to the branes at y — 0, 1 , respectively. We impose Z2 
symmetry across the two branes. That is, for any given function Q, we assume that 

[Q'] = 2Q'(0 + ) , [Q\ = -2Q'(r) . (8) 

To conclude, we describe the four dimensional induced metrics of the two branes. Since 
they are at fixed positions, their induced metrics are simply given by 

ds 2 = -dr 2 + a 2 (r) rfx 2 . (9) 

That is, we recover a FRW Universe with proper time dr = e Bi dt , and scale factor a = e Ai 
(where A4 and B t refer to quantities evaluated at the positions of the two branes). The 
Hubble parameters on the two branes are thus given by 



_ 1 da 
tii= - — 
a dr 



e~ Bi Ai , (10) 



where, as usual, dot denotes derivative with respect to the bulk time t . The Hubble param- 
eters Hi are invariant under residual gauge transformations of the metric 



3 Numerical Code 

In this section, we describe the algorithm that we employ to integrate the equations of 
motion (jHJ) numerically. The algorithm copes with two tasks: It provides the time evolution 
of iV + 1 grid sites, equally spaced between the two branes at y = 0, 1 , and it solves the 
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constraints arising from the boundary conditions at the two branes. In both second 
order discretization scheme is used. In the current version of the program, the same step size 
of discretization is employed in both the time and spatial directions: dt — dy = 1/JV = e . 
This assumption is made not only for simplicity, but also to assure the proper propagation 
of the numerical data along the characteristics of the partial differential equations. 

It is convenient to scale the factor \ = Mf out of the action ffl. This fixes the units 
of the scalar field and its potentials. However, the dynamical equations of motion admit 
a scaling of the metric functions, which allows to choose, in principle, arbitrary units of the 
space-time scales. We use this freedom to secure the supergravity limit of our models, this 
is to say that all length scales are greater than the 5d fundamental scales £5 = M§ . We 
discuss this issue at a greater length in Section 13.41 

3.1 Bulk Evolution 

The system of bulk equations consists of the three second order differential equations (0) 
for the functions A, B, and <ft , which we call bulk evolution equations, as well as the two 
constraint equations (jUJ). The latter are preserved by the evolution equations, and can be 
used as a check of accuracy of the numerical integration. 

In the evolution equations derivatives of functions only appear in the forms / — /" 
and / g - f'g' . 

We discretize the equations by finite-differencing these combinations using the leapfrog 
scheme (see Figure HJ). Let denote the value of the field / at a given grid point y^ on the 
last time step t that had been computed, / nr = / (t, yi) • Then define f\ t and f rt to be the 
value of / at the same time t and on the left and right neighboring sites, f\ t = f (t, and 
/ rt = / (t, yi+i) ■ Finally, denote with / dn and / up the value of the function on y; L at the two 
times just before and after t , respectively /dn = / (t — dt, yi) and / up = / (t + dt, yi) (see 
Figure [TJ. In terms of these quantities, the relevant differential operators of / at the point 
(t, yi) can be discretized with second order accuracy as 

/-/" = \ (/up + /dn - /it - /rt) + O (e 2 ) , 
/ 9 - f'g' = ^ [(/up - /dn) G?up - Sdn) - (frt ~ /it) (fe " 9lt)) + O (e 2 ) . (11) 

Recall that e — 1/N corresponds to the distance between consecutive grid sites. After the 
discretization, the three evolution equations become three algebraic equations, which can be 
solved for the unknown quantities A up , B up , and (f) up . This procedure is repeated at each 
bulk site, leading to the bulk values of the three functions at t + dt , which are then used in 
the subsequent time steps. 
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Figure 1: Numerical evolution scheme. 



3.2 Boundary Conditions 

The numerical scheme described in the previous subsection allows us to determine the value 
of the metric coefficients and of the scalar field at the next time-step for all the bulk sites, 
but not for the two sites i = 0, N , corresponding to the positions of the two branes. To 
obtain the latter, the boundary conditions (J7J) have to be used. First we advance all the bulk 
sites as described in the previous section. Once we know the value of A, B, and i n the 
bulk at time t + dt, equations (jZJ can be finite-differenced into a set of algebraic equations 
for the boundary values at that time. In the following we describe how to implement this 
procedure at y — . The computation for the other brane proceeds analogously. 

The boundary conditions contain first derivatives with respect to y of the metric coeffi- 
cients and of the scalar field at the brane locations. An asymmetric discretization for the 
first derivative of a generic function / , which preserves second order accuracy in e , is given 
by 



fo = ^(-3/0 + 4/x 



h) + O (e 2 ) 



(12) 



Since the right hand side of the first two boundary conditions coincide, we replace the 
first of them simply by A' Q — B' Q = , or, using equation (fT^jl. 



-3 (Aq - B ) + 4 (At - B x ) - [A 2 - B 2 ) = 



(13) 



If we define j3 = e B , the second boundary condition simplifies to (1//3)' = Uq (0o) /6 , which 
can be rewritten as 

a = ZPifo (u] 

Finally, the third boundary condition gives 

4 0i -0 2 -3 0o -eft Wo) = • (15) 
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Only the values of the three functions on the brane are unknown. By substituting the 
value for /3 given by equation (|14|) into equation (JT5J), the latter becomes an equation 
where the only unknown quantity is 0o • For specific brane potentials Uq , this equation can 
be solved analytically; more generally, one can solve it numerically through some iterative 
method. In our algorithm, the iterative Newton's method is employed. Finally the value of 
B = ln(/3 ) can be used in (fT^j) to determine A . 

3.3 Initial Configurations 

Initial conditions are imposed by specifying the three functions and their first time derivatives 
on the grid sites at some initial time t = . We denote them as A . . . , O (jji) > with 
i ranging from 1 to N — 1 in the bulk (i = 0,N are the sites of the two branes). These 
functions can not be chosen arbitrarily but rather must satisfy the constraint equations (JHJ). 
Once this is done a second order Runge-Kutta time step is used to "convert" the initial 
conditions of the form f (y), f (y), into initial conditions given at the first two initial time 
steps fo (y) and fo+dt (y)- The Runge-Kutta step is done as follows 

fo+dt (y) = fo (y) + dt (j (y) + X -dt f (y)j , (16) 

where the second time derivatives fo(y) are replaced by the equations of motion (JHJ). This 
"conversion" is needed to have the initial conditions in a form suitable for the leapfrog scheme 
described in Subsection 13.11 

In general, the initial time derivatives can be non-vanishing, so that one can study situa- 
tions in which the geometry of the extra dimension is time-dependent already at the begin- 
ning of the numerical integration. For example, this is the case for the AdS-Schwarzschild 
solution we will deal with in Section 14.21 For each such case the choice of initial conditions 
must be consistent with the constraint equations (|Sj). We discuss one such algorithm in 
Section IP 

A particularly interesting class of initial conditions is however the one of static warped 
solutions 

ds 2 = W(y) 2 (dy 2 - dt 2 + e 2Ht dyi 2 ) . (17) 

characterized by a fixed bulk geometry and maximally symmetric (de Sitter or Minkowski) 
branes. This metric turns to the form (jSJ) with the identification 

B (t, y) - B (y) = In W , A(t,y) - B (y) + Ht , (18) 

where H is the Hubble parameter of the de Sitter brane and the bulk scalar field (f> is also 
a function of y only. Such solutions were studied with dynamical system methods in |17j . 
The numerical integration can be used to check their stability. Numerical errors due to 
the grid discretization act as small perturbations. If the initial configuration is not stable, 
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the tiny numerical errors accumulate with time, and eventually lead to an evolution of the 
system. When this happens, a full numerical calculation is the only tool to study where 
this evolution leads to, namely whether the two branes collide, move apart to infinity, or 
get stabilized at a finite distance in another static but stable configuration. As we will 
see below, in many cases de Sitter branes turn out to be unstable. Therefore even static 
warped geometry configurations can provide suitable initial conditions for time-dependent 
braneworld dynamics! 

When numerical inaccuracy is used to seed the evolution, as described above, the initial 
amplitude and consequently the timing of the instability depends on the accuracy of the 
numerical integrator. This accuracy is in turn related to the spacing of the grid sites in 
the bulk. Increasing the number of grid sites decreases their separation, and the instability 
develops later. Alternatively, initial perturbations on the top of the static configurations 
can be imposed directly as initial conditions. This allows a quicker development of the 
instability, or, for static configurations, the excitation of some of the lowest eigenmodes 
of the system. A simple class of initial perturbations, which can be implemented in our 
numerical algorithm, is described in Appendix IB. II and illustrated in Figure 03 We found, 
however, that the qualitative behavior of the system did not depend on the details of how 
the initial perturbations are generated, whether imposed explicitly or through numerical 
roundoff errors. In the following, we therefore discuss instead how the static configurations 
are determined. 

For static configurations, the bulk equations reduce to 

0" + 3 B' 0' - e 2B V (0) = , 

B ' 2 v 2By ^-^' 2=ij2 ' (19) 

In addition, the last two of the boundary conditions (J7J) have to be satisfied at each brane. 
In the gauge we are using here the constraint equations © are automatically satisfied. 

The bulk equations are thus reduced to a system of first order differential equations for 
the functions -8,0, and 0' , so that the phase space of possible solutions is effectively three- 
dimensional [17J. To solve these equations subject to (given) boundary conditions, we specify 
the values of the three functions at y = , as well as the value of the constant parameter H*. 
For a given brane potential Uq , only two of these four numbers can be chosen arbitrarily, and 
the other two are determined by the junction conditions at the first brane. (In the 3d phase 
space this means that the junction condition at the first brane defines a Id curve in phase 
space along which the trajectory must begin.) The bulk equations (j!9|) are then integrated 
with a standard fourth order Runge-Kutta integrator. Depending on the initial values and 
on the bulk potential, the bulk solution may become singular before the brane at y = 1 is 

*One may be wondering why we can specify four variables in a 3D phase space. Recall, however, that 
in our gauge the position of the second brane is fixed at y — 1. In the language of ^7] we are using three 
degrees of freedom to specify the starting point of our trajectory in phase space and one to specify the length 
of the trajectory, i.e. at what point on the trajectory the second brane will be found. 
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encountered. If this happens, some other initial values (or some other bulk potential) have 
to be considered. 

Even if the brane at y — 1 is reached, we face the nontrivial problem of satisfying the 
boundary conditions also at the second brane. The simplest way to solve it is to regard 
the junction conditions as equations for the parameters of the brane potentials. One can 
freely choose the three numerical values at the first brane (as well as the numerical value 
of H), integrate the bulk equations, and then use the junction conditions to determine the 
potentials at the two branes. ' However, one is typically interested in the more difficult 
situations in which the brane potentials are specified, and the initial configurations have to 
be determined accordingly. 

In the second case, we face a boundary- value problem: values of the fields satisfying the 
boundary conditions at the first brane do not in general lead to field values that satisfy the 
boundary conditions at the second brane, once they are evolved across the bulk according the 
bulk differential equations ()19|) . It is by no means guaranteed that any choices consistent with 
the junction conditions at both branes exist. Indeed, as discussed in [T7j, many potentials 
do not give static solutions at all, while some others typically lead to only a finite number of 
them. In Appendix El we discuss the numerical method (known as the "shooting" method 
[22]) which we employ to find these solutions. 

3.4 Units 

Let us inspect the dynamical equations ©, constraint equations ©, and the boundary 
conditions (J7J). While the units of the bulk scalar field are fixed by our form of the action (fT|). 
it is easy to see that these equations are invariant under the following scaling transformation 

A^A + S', B ^B + S, (20) 

where S' and S are arbitrary real valued transformation parameters. The scalar field poten- 
tials enter the equations only in the combinations e 2B V and e B U. Therefore (j2T)j) can be 
accompanied by the transformation 

V^e~ 2S V, U^e~ s U. (21) 

Suppose some metric functions A, B are the solutions of © for given potentials V, U. 
The scaling transformations (j2*Uj) and (}2*T]) tell us that from these metric functions we can 
generate a family of solutions for re-scaled potentials. This is very useful for introducing 
the units of scales for numerics. Indeed, while y and t in (J3J) are 2d conformal length and 

tin general, this does not determine the brane potentials, but only their values and their derivatives at 
a single value of the field <j> . One can complete the functional form of the potential arbitrarily, say as in 
equation J3J) • In this case, one is for example free to choose large positive values for the two mass parameters 
Mi , favoring values of the scalar field at the branes which are close to the vacuum expectation values Oi . 
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time (i.e. affine parameters along corresponding directions), the metric function e B defines 
the physical interbrane distance D and the physical time. As often occurs in numerical 
simulations, it is not always easy to extend the range of variables, like e B in our case. As we 
will see in the example of the next Section, numerical stability (without brane stabilization) 
has a controlled but finite life time. If we naively increase the scale of e B , the stability will 
be short-lived. The trick is to continue to work with numerically convenient values of e B , but 
interpret scales in units of / = e s I5. One can take S to be large enough to have the scale e s I5 
much greater than the fundamental bulk scale I5. This is to say that numerically we solve 
our equations not only for a given scale and given choice of parameters of the potential, but 
for the whole family of scales and parameters which corresponds to the orbit of the group 
transformation (|20|) . (j21j) . For the parameters of the potentials (J2J) we have the following 
units: [m] = e~ s M 5 , [M] = e" 5 M 5 , [A] = e^Mf , [A] = e~ s M^ [0] = [a] = M~ 3/2 . 

The time evolution of variables in the paper will be plotted versus conformal time t. 
The units of t are the light crossing time between branes. This corresponds to the distance 
between branes in the conformal coordinate y which is simply 1 in our units. 

As usual, the parameters for the numerical simulations do not allow the introduction of 
a large hierarchy, since numerical inaccuracies accumulate much faster. Therefore most of 
the parameters are chosen to be of order unity in our units. The values of parameters that 
correspond to the numerical runs shown in the figures of this paper are collected in the Table 
[I] in the Appendix. 



4 Branes in AdS Background Without a Scalar Field 

The algorithm presented above allows an exact integration of general two-brane configura- 
tions with bulk scalar fields. In this and in the next sections we discuss in detail several 
applications. The first two examples of this section have no scalar field, and the evolution 
is known analytically. We report them mainly to discuss the accuracy of the code and to 
outline its main features. The code is accurate enough to reproduce the known analytic so- 
lutions. In the following sections we will study the more complicated non-linear evolution of 
a system with a scalar field, for which the solutions were not previously known. Fortunately, 
we still will be able to check certain properties of the solutions analytically. 

In this section we first consider the static unstabilized Randall-Sundrum flat brane so- 
lution, from the point of view of the numerical solution of the equations. Then we study 
non-static (moving) branes in an AdS-Schwarzschild background, and compare the numerical 
solution with the known AdS-Schwarzschild solution. 
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4.1 The Randall-Sundrum Model 



Our first example is the two-brane Randall-Sundrum model It represents a particularly 
simple example of a brane world that only consists of a five dimensional AdS space with a 
curvature radius I 2 = — 6/A, determined by its 5d cosmological constant A, and of two flat 
branes with tensions Aj = ±6/Z . The system is entirely described by one time independent 
function. In terms of the 2d conformal gauge Q we have 



A{y)=B{y) = -\n 



y+(e D / l -iy 



(22) 



where D can take any constant value, which - according to equation (0J) - corresponds to 
the interbrane distance. We can reproduce this set-up in our code by simply setting to zero 
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Left Panel: Drift of the free modulus of the RS model depending on the numerical 
Right Panel: Stabilization of the RS radion. 



all the scalar field related parameters in the bulk/brane potentials in equations (jSJ) (as well 
as the initial conditions for the scalar field). The numerical solutions of the equations (0) 
are in agreement with (|2*2*|) . As discussed in section EHfl small perturbations are unavoidably 
introduced by the discretization. Thus, this set-up is particularly useful for verifying the 
accuracy of the code. Notice also, that numerical instability is much worse for the RS 
without stabilization. In the left panel of Figure El we show how the time-scale at which the 
instability develops is related to the number N of bulk sites. The more we increase N the 
more the accuracy of the computation increases, and numerical instability is delayed for the 
later times. We estimate the time scale where the code is stable as being proportional to 
the grid resolution N . 

The right panel of Figure |2] also shows how the introduction of the stabilization mecha- 
nism (with the bulk scalar fields with the potential (J2J)) can, for appropriate choices of the 
parameters, lead to a stabilization of the interbrane distance. In this his case the code is 
much more stable. We discuss this issue in more detail in Section O 
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The choice of parameters and initial conditions that was used in the numerical runs 
plotted in Figure 121 as well as the ones for all following simulations are collected in Table ^ 



4.2 The AdS-Schwarzschild Solution 

Starting from a setting similar to the Randall- Sundrum example of the previous section, but 
allowing for non-vanishing initial (at t = to) time derivatives, we generate time-dependent 
numerical solutions that belong to a larger class of solutions. Assuming the initial spatial 
profile A (t Q ,y) = B (t , y) of equation (J2~2~j) . the constraint equations © are solved by 



A(t ) = c 



1 

y + 



D/l _ 1 



B(t ) = -A{t ) , (23) 



where c is a constant. The choice c = gives Randall- Sundrum solutions, while a non- 
vanishing c corresponds to moving branes. From the Birkhoff theorem for plane-parallel 
brane configurations it follows JH] that the generic 5d bulk metric must be a stripe of 
the AdS-Schwarzschild geometry (where the Schwarzschild mass is virtual because it is 
screened by the branes). Thus, the branes are moving in an AdS-Schwarzschild background. 
To see this, note that in the absence of the scalar field and for brane tensions as in the 
Randall-Sundrum model, the boundary conditions (J7|) give A' (A' + B') = 2Aexp(2£>) at 
the location of the two branes. From the bulk equations, we then recover 

A + 2A 2 -AB = , y = 0,l, (24) 

which, in terms on the proper time r and the Hubble parameter Hi on the two branes (cf. 
equations © and (JTUJ) ) , becomes 

^. + 2Hf = 0. (25) 
This corresponds to a radiation dominated standard four dimensional Universe, 

The appearance of effective radiation domination on the branes is characteristic of an AdS- 
Schwarzschild bulk geometry ^H]- The invariant of the 5d Weyl tensor C 2 = CabcdC abcd 
projected into the brane scales as C 2 oc a~ 8 , where a (r) is the scale factor of the induced 
FRW brane metric. Since a (r) is radiation dominated, at the brane we have 

C 2 = C 2 [l + 2H (r-T )]- A . (27) 
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Dark Radiation Weyl Tensor C =C ABCD C 
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Figure 3: Comparison between numeric and analytic solutions for the Hubble parameter and 
Weyl tensor on the brane. 

The two equations (J2fijl and (|27j) are independent of the choice of coordinates in the bulk, 
and can be easily reproduced with our code. For a particular realization of the initial condi- 
tions ((SHI with the parameter c = 1, in Figure El we plot the numerical calculation (squares) 
of the time evolution of H and the 5d Weyl tensor on the brane. Solid curves correspond to 
the AdS-Schwarzschild analytic solutions (ffijj) and (j27jl . The agreement between numerics 
and analytics is manifest. 
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Figure 4: Left panel: Non-physical waves appearing as gauge modes in the metric function 
B (t, y) (for c = 1). Right panel: interbrane distance for various initial conditions controlled 
by the parameter c. 



In the left panel of Figure EJ we show instead the evolution of the metric component B (t, y) 
for the same configuration used to generate Figure El The ripples of B(t, y) are not physical. 
As mentioned in Section |21 our choice of coordinates does not fix the gauge completely. The 
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residual gauge freedom appears numerically as ripples in B(t, y). The precise form of these 
gauges modes is worked out in Appendix IA.21 and is in agreement with the numerical plots. 
The lowest frequency mode of these gauge modes generically appears in the evolution of 
B (t, y) . In the left panel of Figure EJ we see this effect in the form of two bulk waves with 
period 2 , which propagate on top of the profile of B (t, y) . As discussed in Section EJ these 
gauge modes do not affect the interbrane distance, as defined in equation (jlj). 

The numerical evolution of the interbrane distance for the AdS - Schwarzschild solution 
is shown in the right panel of Figure 0] for different values and signs of the parameter c. For 
positive c the branes approach each other, while for negative c they move apart. 

5 Instability of de Sitter Branes and Restructuring of 
Warped Configurations 

Let us now study the evolution of the system J5J) in the presence of a bulk scalar field. We 
have to specify initial conditions, which do satisfy the constraint equations. This task is now 
more complicated than it was without the scalar field. 

5.1 Instability of Warped Geometry with Curved Branes 

As a starting point we can check the code for known static warped geometry configurations 
()18|) with the scalar field and potentials chosen so as to stabilize the branes. We can then 
impose perturbations consistent with the constraint equations. (This technique is described 
in detail in the Appendix IB.ljl 

Static solutions of warped geometries with bulk scalar fields and with branes at the 
boundaries have been studied and classified in In the 2d conformal gauge the static 

solutions with curved branes are given by 

A(y,t)=B(y)+Ht ,4> = 4>(y), (28) 

where B(y) and (f>(y) are related through a set of ordinary differential equations, which can 
be treated with the methods of [T7j . 

We use scalar field potentials (|5Jl. which are designed for brane stabilization. The outputs 
of numerical integration of an initially static configuration of two curved (de Sitter) branes 
and bulk scalar filed with small perturbations around it is shown in Figure EJ We see the 
appearance of time-dependence in the initially static field 4>(y), departure of metric function 
A from the hypersurface described by equation A(y,t) = B(y) + Ht, and a change in 
the interbrane distance D. We show two realizations of this model with different initial 
perturbations. From these results we conclude that, surprisingly, the static solutions with 
scalar field potentials that are supposed to stabilize branes are unstable for a range of H\ 
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4>(t,y) 






Figure 5: Instability of static solutions with dS branes: Perturbations have induced sig- 
nificant departure from the static solution at t ~ 40. The two unstable solutions shown 
correspond to positive (blue, upper surfaces and increasing D) and negative (red, lower 
surfaces and collapsing D) initial perturbations of 5<ft, see Appendix lB.il 



Fortunately, this unexpected result, which we found here numerically, can be indepen- 
dently obtained with analytical methods reported in the accompanying paper [THj. Indeed, 
it is possible to consider linearized perturbations of the bulk scalar field 5(f) and scalar metric 
perturbations 

ds 2 = W(y) 2 [(1 + 2®)dy 2 + (1 + 2*) {-dt 2 + e 2Ht d* 2 )] . (29) 

around the background warped geometry (|17|) . where $ and $ are small metric perturbations. 
From the linearized Einstein equations one can derive second order differential equations for 
the fluctuations, which can be factorized into 4d massive scalar harmonics on the de Sitter 
slices and KK eigenfunctions with eigenvalues m. The lowest eigenvalue in the KK spectrum 
corresponds to the radion mass. The lowest eigenvalue m 2 is estimated as 

m 2 = -AH 2 + m 2 (H) , (30) 

where mo = | j ;j~^7=2 is a functional of H. In many cases the first (negative) term in 
()43j) exceeds the second positive term, causing a tachyonic instability of the curved branes. 
Indeed, the temporal part of the eigenmodes has an exponential instability 

f m (t) oc e"', (31) 

where 

/ [9 im^T 3\ TT 

"=(V4 + W-2j ff - < 32 > 

The time dependence from numerical calculations corresponding to Figure El is consistent 
with the analytic result (}3Tf . In the limit of H = formula (J30|) is reduced to the known result 
for flat branes where the branes configuration is stable ^H]- The curvature of the branes 
upsets the balance between the bulk scalar gradients and its potentials, which otherwise 
provide stabilization. 
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Thus, both from numerics and analytics we conclude that many static configurations 
with de Sitter branes are unstable against classical (or quantum) fluctuations. While in the 
following we mostly discuss the physical meaning and consequences of this result, here we 
also note that this effect provides us with a tool to study brane dynamics numerically. This 
is because we can start with a simple controllable static configuration, without needing to 
resolve the time- dependent constraint equations. Doing so, we can investigate numerically 
fully non-linear time dependent dynamics due to the real physical tachyonic instability of 
the initial configuration. 

Depending on the sign of the initial perturbations (the coefficient p in equation 
we encounter runaway behavior towards smaller or larger interbrane distances as shown in 
Figure 03 We consider this type of non-linear dynamics in Section |3 Sometimes we do not 
find a runaway behavior, but rather a re-structuring of brane configurations as a transition 
between (at least) two static warped geometries. This case will be considered in the next 
subsection. 

5.2 Dynamical Transition between Two Static Solutions 

As we discussed in the Introduction, the construction of brane models with de Sitter branes 
is particularly challenging. Stable static solutions with inflating branes can only be achieved 
provided the spatial gradient of the bulk scalar field is sufficiently high, cf. equation (J30J). 

In the context of static warped geometries, brane embeddings can be investigated in 
geometrical terms in a three dimensional phase space J7| . This technique is especially useful 
to show that more than one static solution for a given brane model, i.e. given potentials V(4>) 
and Ui(4>), might exist as illustrated in Figure El Many of these solutions are unstable, as 
shown above. A fully numerical integration is a powerful (and maybe the only) tool to study 
the non-linear dynamics of the unstable brane configurations. A more comprehensive study 
of this issue, with different bulk/brane potentials taken into account, will be presented 
elsewhere. Here we limit our discussion to potentials of the class (0). In this subsection 
we discuss the case in which the braneworld model admits one unstable and one stable 
static solution, and the evolution of the system drives a transition between the two. Small 
perturbations around the unstable solution trigger the tachyonic instability of the system, 
which is followed by a rapid evolution of the bulk configuration. 

An example of dynamical transition between two static brane configurations is shown 
in Figure |H1 We plot the time evolution of the distance D(t) between the branes, the 
radion mass m 2 , the Hubble parameter (curvature) H, and the Weyl tensor invariant C 2 . 
The last two are defined as the averages of these quantities over the extra dimension. We 
observe a transition between two states, from an initial brane configuration with higher 
brane curvature (larger H) to a final configuration with lower curvature (smaller H). The 
first state is unstable; during this regime the radion mass is tachyonic, m 2 < 0. The value 
H decreases with time until the second term in ()43|) dominates and the tachyonic instability 
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ceases. In the cases we have studied, the decrease of H is accompanied by a decrease of 
the physical interbrane distance, until the stable configuration is reached. * The final static 
configuration has positive m 2 . 

The dynamics of the transition between the two static configurations is quite violent, and 
is accompanied by a burst of the Weyl tensor C 2 . The value of C 2 vanishes for the warped 
geometry configurations at the beginning and the end of the transition. 
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Figure 6: Left panel: Transition between two static brane configurations. Right panel: 
Transition observed on the brane at y = 1 for various values of the parameter M of equation 
(J2J), where the last plot M = 1 corresponds to colliding branes (see Section EJ). 



Remember that we restrict ourselves to (t, y) dependence and "planar" symmetry of the met- 
ric. Of course, the actual dynamics between two warped configurations does not necessarily 
occur in this class of metrics, and 3d inhomogeneities along the brane can be excited. As 
shown in the tachyonic instability of warped geometry with de Sitter branes occurs for 
scalar long-wavelength inhomogeneous modes with 3d momenta k. The present form of the 
BraneCode cannot take them into account. We assume that the background k — > mode 
dominates, but this should be investigated in the future. Tensor inhomogeneous modes do 
not have tachyonic KK spectra [21] around the curved brane warped geometry. In fact, 
gravitational waves are absent for systems with planar symmetry. However, based on the 
evolution of the Weyl tensor C 2 , we conjecture that the actual dynamics is accompanied by 
a burst of gravitational wave emission with 3d momenta of the order of the non-adiabatic 
frequency ~ 1/At, where At is the time of transition. It will be interesting to check this 
with a linear tensor mode analysis around the background geometry of the Figure El 

■'"The quantity H was defined in equation l|18(l only for static configurations. During the time evolution, 
we choose to dehne it as the average over y of B (t, y) — A (t, y) at any fixed time t . In the examples 
discussed, we saw that the combination A (t, y) — B (t, y) depended only weakly on y during the whole 
evolution. 
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It is interesting to follow how the final state of the unstable warped configuration depends 
on the parameters of the potentials. We illustrate this with the parameter M of the potential 
(J2J). In the example shown, the four dimensional cosmological evolution on the two branes is 
characterized by a transition between two de Sitter spaces. In the right panel of Figure El we 
show how the four dimensional Hubble parameter on one of the two branes changes as we 
change the brane mass parameter M of the scalar field, while leaving the other parameters 
unchanged. In the limit of large M, the value of the scalar field on the branes is always 
very close to its expectation value a . Moreover, phase space portraits [T^j (see Appendix lB|) 
indicate that the two static configurations approach each other in this limit. This is also 
visible in Figure® where we see that the difference between the initial and the final value of 
H decreases as M is increased. As an analogy, one may say that higher mass parameters M 
correspond to more rigid systems, characterized by stiffer and quicker transitions between 
the two static regimes. Tuning the parameters of the model, one can have flat Minkowski 
branes in the stable final configuration. 

In the limit of negligible M the system does not admit stable configurations at all. The 
curve with M — 1 , shown in Figure |H1 corresponds to a case in which the dynamics of the 
system lead to a collision between the two branes. This case is discussed in detail in the 
next section. 

6 Brane Collisions 

Unstable warped configurations of curved (de Sitter) branes provide suitable initial condi- 
tions for studying colliding branes, as we saw in Subsection 15.11 By controlling the initial 
fluctuations (see Figure |SJ) we can generate numerical runs with brane collisions. 

The collision of branes is an interesting subject by itself. In cosmology colliding branes 
appear in models of brane inflation [6 as well as in models without (early universe) inflation 
[23] . The latter models have difficulties which were discussed elsewhere (sec e.g. [26 ). In 
this paper, we focus on the issue of the bulk geometry and scalar field profiles of colliding 
brane configurations, in a more general context. 

In the next Subsection I6.1l we show a numerical example of the brane collision and try to 
understand the properties of the interbrane geometry. We find that they become independent 
of the specific brane/bulk potentials of the model. In Subsection 16.21 we further argue that 
there is a universal Kasner-like space-time asymptotic of the interbrane geometry. This is a 
strong-gravity regime which can not be described in 4d by the moduli approximation. 

6.1 Geometry Between Branes 

Figure d and Figure |H1 show in detail the evolution of the bulk scalar field (f)(t,y), metric 
functions A(t,y), B(t,y) and interbrane distance D(t) in runs which begin with an unstable 
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warped de Sitter brane configuration and end with a brane collision. The first thing to 




Figure 7: Flattening of 4> gradients during the brane collision. The left panel shows 4>(y) for 
different t, going from top to bottom. The lower panel shows (f>(t) different y. In either plot 
you can see that <fi becomes nearly homogeneous at later times. 

notice is that the system becomes homogeneous along the y coordinate. This is seen as the 
flattening of gradients over time. A similar flattening in y direction occurs for the metric 
components, see Figure |H1 Also notice that the absolute value of <fi increases with time. This 
increase can be fit well by 4>(t) ~ Int. 

A second feature of the brane collision is the decrease of the metric component e B ; 
asymptotically B — > — oo during the collision, cf. the definition (0J) of the interbrane distance. 
Recall that the bulk/brane scalar field potentials in the bulk equations © and boundary 
conditions (JJJ) are always multiplied by exponents e B . Therefore the contribution of the 
bulk/brane potentials becomes more and more negligible in the dynamical equations (0) 
during the collision. 

This leads us to the important conclusion that asymptotically the dynamics of the brane 
collision do not depend on the form of the bulk/brane potentials. Notice, however, a poten- 
tial exclusion from this rule related to exponential potentials e°"^. In this case the typical 
logarithmic time divergence of <fr leading up to the collision leads to the growth of the value 
of the potentials with time which may compensate the decrease of the metric function e B . 
In this paper we concentrate on the potentials (J2J) where asymptotically the dynamics are 
potential-free. 

This potential-free asymptotic immediately helps to explain heuristically the first feature, 
why the system becomes homogeneous along the y coordinate. Indeed, looking at the bound- 
ary conditions, we see that the gradients of A , B , and at the branes are proportional to 
exp(B) and therefore vanish as e B — > 0. 
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Next, let us consider equations (jHJ) under the assumption that the bulk/brane potentials 
for the scalar field can be neglected and that the geometry becomes homogeneous. After 
this simplification equations (j5J) become ordinary differential equations, which can be easily 
solved. We find 

A = A + -\n(t c -t) , 
cj) = O - r In (t c - 1) , 

B = fi - ^ (l - 3i?ot c - f r 2 ) ^ - i (1 - |r 2 ) ln(t c -t) . (34) 

The constants of integration A ,B Q , O correspond to the values of the fields at some time 
t = 0. The time t = cannot be the beginning of integration where we know that the 
approximation does not hold. We will give meaning to the integration constants shortly. 
The collision time is t c = —1/Aq. We also introduce a convenient intermediate parameter 
t — (pot c - The brane collision corresponds to values of r satisfying r 2 < |. The scalar 
field potentials, as well as inhomogeneities along the y coordinate, result in small corrections 
which are neglected in ()34|) . Asymptotically the logarithmic terms in (J34|) dominate and we 
arrive at a homogeneous metric with power-law dependence on time. This is nothing but 
the recognizable Kasner-like space-time metric. 



6.2 Universal Kasner-like Asymptotic 

The regime when the logarithmic terms in ([3*4*]) determine the behavior of the system corre- 
sponds to the universal Kasner solution in five dimensions with a massless scalar field. Four 
dimensional homogeneous but anisotropic Kasner solutions with the massless scalar field were 
constructed a long time ago in [27] . Its higher dimensional generalization is obvious [28j. 
Indeed, in 5d we have the following exact solution with the massless scalar field 

3 

ds 2 = -dr 2 + r 2p ydy 2 + r 2pi dx 2 , 

i=i 

Pi + P2 + P3 + Py = 1 

p\ + pI + pI+pI = i-9 2 

= glnr . (35) 

The vacuum Kasner solution has q = 0. The parameter q characterizes the contribution of 
the scalar field. The time t in the 2d conformal gauge © and r are related by transformation 

t = T l ~ p y. (36) 

The significance of the Kasner-like space-time (135)1 is not only in the fact that it is an exact 
solution of the Einstein equations, but mostly because it is a generic asymptotic of arbitrary 
collapsing solutions [25] . 
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In this section we explicitly demonstrate how the geometry of colliding branes, as a case 
of the collapsing solution, approach the universal Kasner-like asymptotic. 

Kasner-like geometry as generic collapsing solution was already advocated in string cos- 
mology [SH]. As we show here, this asymptotic also applies to brane cosmology (in other 
words string cosmology with branes). There is, however, a specific new feature that ap- 
pears in the brane cosmology case. The isometry in the brane directions is reflected in the 
additional constraint. In 5d 

Pl=P2=P3 ■ (37) 

This constraint and the two equalities for Kasner indices allow to express p\ and p y through 
the parameter q 

pi=i(i±\/i-F) ■ p * = \( i * 3 f-P)- <38) 

The range of the parameter gis— ^<g<^; the range of p y and p\ correspondingly are 
~ \ < Py < 1 an d < pi < 7j. In the vacuum limit of vanishing q one finds p\ = or 1/2 
and p y = 1 or —1/2. 

One of the feature of the Kasner-like asymptotic is a chaotic alteration of the indices 
Pa |29j . with alternating contraction and expansion in some of the directions. In the presence 
of the scalar field the process ceases as all pa become positive. For colliding branes p y > 0, 
which gives < p\ < |, so both indices are positive and no alteration of indices is expected. 

For our case (jUTJ) the Kasner solution can be rewritten to the 2d conformal gauge 
with the help of the time redefinition (J3l)j) 

ds 2 = (dy 2 - dt 2 ) + tUx 2 . (39) 
In terms of the metric functions A and B and the field 0, the Kasner solution (|3*9*j) reads as 

A = l -\n{t c -t) , 

B = -^-hx{t c -t) , (40) 

J- Py 

= ln(t c - t) . 

l-Py 

The solution (J4*U|) is identical to the leading terms of ()34|) by the identification q = 6r/(3r 2 +4) 
Thus the integration constant r in (|34|) is related to the parameters of the Kasner solution. 
Figure |S1 shows how the metric components A and B and the field found numerically, 
approaches the universal Kasner asymptotic (}3T?|) . Next, consider metrics induced by the 
bulk Kasner geometry at the branes (which is independent of y) 

ds 2 = -dr 2 + (r c - r) 2pi dx 2 , (41) 
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Figure 8: Numerical solutions (lower red surfaces) asymptotically approach universal Kasner- 
like solution (upper blue y - independent surfaces). 



The induced Hubble parameter on either branes is then given by 

Pi 



H 



t c -t 



(42) 



This time-dependence of the Hubble parameter at the brane is a good fit to the asymptotic 
behavior of the Hubble parameter we found numerically, as illustrated in Figure EJ The 
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Figure 9: Induced Hubble parameter H(t) at the colliding branes 

induced metric at the brane (|41j) depends on the parameter q through the index p\. This 
parameter is absent in the simple moduli approximation of the 4d effective theory at the 
brane, which does not take into account strong gravity arising in the bulk geometry. 
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7 Summary 



We design the numerical code, BraneCode®, to treat non-linear time-dependent dynamics 
of 5d braneworlds with plane-parallel branes at the edge and with bulk scalar field with 
arbitrary bulk and brane potentials. It is possible to choose a convenient gauge where 
the brane positions are not changing with time, and dynamics is imprinted in two metric 
components and bulk scalar field. These bulk equations for gravity and the scalar field are 
supplemented by boundary conditions at the orbifold branes and initial conditions in the 
bulk. We also treat the constraint equations at the initial time hypersurface. So far we 
have only included a single bulk scalar field, but the code could in principle be extended to 
include other layers such as additional scalars in the bulk or on the branes. 

We check the code for the brane models with known analytic solutions. For two branes 
embedded in the 5d background with negative 5d cosmological constant without the scalar 
field we numerically reproduce generic AdS-Schwarzschild solution 

Next, we consider more comprehensive braneworlds with the bulk scalar field. We in- 
vestigate numerically small perturbations around warped stationary configuration with bulk 
scalar and with de Sitter branes including the bulk/brane potentials, which are introduced 
for the brane stabilization. However, for the large enough 4d curvature of inflating branes 
the system is unstable and runs away from the initial warped configuration with de Sitter 
branes. This effect is confirmed independently by an analytic calculation of small scalar 
perturbations in this setting [T^j. The scalar fluctuations around a warped configuration 
with curved branes have as their lowest eigenvalue 

m 2 = -AH 2 + m 2 (H) . (43) 

The term m (H) is a functional of H, and depends on the parameters of the model. If 
parameters are such that m 2 becomes negative due to excessive curvature ~ H 2 , the brane 
configuration becomes unstable. For relatively low values of H 2 the radion mass ()43j) is 
positive and the system is stable. Our interpretation of this instability is the following. 
Stabilization of flat branes is based on the balance between the gradient <$' of the bulk scalar 
field and the brane potentials U(<p) which tends to keep <fi pinned down to its values <pi at 
the branes. The interplay between different forces becomes more delicate if the branes are 
curved, and for the brane curvature exceeding some critical value the brane configuration 
becomes unstable. 

Tachyonic instability of curved branes has serious implications for theory of inflation in 
the braneworlds and is discussed in details in the accompanying paper [TKj . 

Our numerical simulations allow us to follow dynamics of the brane configuration trig- 
gered by the tachyonic instability. The end point of evolution depends on the presence or 
absence of another (or others) warped stationary configurations in the model. 

The question about the multiplicity of warped solutions can be studied in the framework 
of warped geometry with no time-dependence. We implemented the geometrical construc- 
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tions in the phase space of solutions of the gravity plus scalar system that had been developed 
earlier. In the model with quadratic bulk and brane potentials, depending on the parameters 
there are single or double warped solutions. 

Thus we see that in some cases the warped branes system can admit two solutions for 
the same parameters of the potentials, with different values of the curvature of the de Sitter 
brane, which is proportional to H 2 . Suppose we start a numerical run with the warped 
solution that has a larger value of the brane curvature, which is unstable. Then we observe 
numerically that this configuration evolves dynamically and ends up in the the state which 
corresponds to the second warped solution. The second solution is stable if the corresponding 
radion mass m? is positive, as in the example shown in the text. This re-structuring is 
accompanied by strong dynamical features like a burst in the Weyl tensor, which vanishes 
in the initial and final warped configurations. Although inhomogeneous tensor modes are 
not included in the code, based on this behaviour of the Weyl tensor we conjecture that 
brane re-structuring should be accompanied by the emission of gravitational waves due to 
the non-adiabaticity of the process. 

All together, this process looks like a decay of the meta-stable state of the strongly 
curved branes due to the tachyonic instability into the more stable state where the branes 
have lower curvature. This transition is marked by a burst of gravitational field anisotropy 
(gravitational waves?). It will be interesting to investigate what applications this may have 
to cosmology with branes. Another potential application of brane restructuring would be 
the problem of the 4d cosmological constant in the braneworld picture. The cosmological 
constant problem was discussed recently from a braneworld perspective, in which a low 4d 
cosmological constant corresponds to a flat brane. There was a suggestion that the flat 
brane is a special solution of the bulk gravity/dilaton system with a single brane [32 j, but 
the model has difficulties [33] • in our setup, we consider two branes. The new element which 
emerges from our paper is the instability of the curved branes. So far we have only shown an 
example of re-structuring between two curved brane configurations. It will be interesting to 
see if there are brane models with more than two stationary warped geometry configurations, 
or with several scalar fields, and to investigate if there is a mechanism for brane flattening. 

Finally, we studied the geometry of colliding branes. As initial conditions we used the 
unstable curved brane configurations with parameters which do not allow another warped 
geometry configuration. In such cases the end point of the brane dynamics is either a brane 
collision or branes moving apart. We investigated in detail the geometry of colliding branes. 
The bulk metric and the bulk scalar field become homogeneous, i.e. y-independent, and 
the brane dynamics asymptotically does not depend on the scalar field potentials. Instead, 
the geometry of colliding branes asymptotically approaches a universal Kasner-like solution 
with a free scalar field. It is known that the Kasner asymptotic is a generic solution of 
the high dimensional gravity/dilaton system [30 J. In our case isometry of the brane slices 
guarantees equality of the Kasner indices p\ = P2 = Pz- 111 5d this condition leaves only 
one single parameter of the Kasner-like solution q, associated with the bulk scalar field 
contribution. This parameter is determined by the initial conditions. For the 5d brane 
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system we considered, there is no chaos in the alteration of the Kasner indices. It will be 
interesting to investigate this issue for other situations, for example when the form field is 
included and the brane dimensions and co-dimensions are different. 
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Appendices 

A Choice of Gauge 

The system we are studying has a gauge freedom which amounts to different possible coordi- 
nates for the five dimensional metric and for the positions of the two branes. Our choice not 
only aims at simplifying the equation of motions, but also at removing "redundant" degrees 
of freedom which would not allow us to write a closed system of equations for the numerical 
integration. In Section we claimed that it is always possible to choose a system of coordi- 
nates in which the (t, y)— part of the metric is conformally flat and in which the two branes 
are fixed at the positions y = and y — 1 , irrespective whether their physical separation is 
constant or changing in time. We show this explicitly in Subsection IA.1I This choice does 
not fix the gauge completely, however. The form of the remaining gauge degrees of freedom, 
which are expected to affect the numerical solutions, is worked out in Subsection IA.21 



A.l Comoving Branes 

By assumption, the system is homogeneous and isotropic along the spatial coordinates x , 
and the position of each brane in the extra space is specified by a function of time only 
(parallel branes). Since the metric coefficients depend only on the two coordinates t and y , 
the metric can be written in the 2d conformal gauge Q 

The change of coordinates 

t f _ f(t + y) + g(t- v ) 

y - , = /(' + »> 2 -'(«-»> , (44, 
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where / and g are two arbitrary scalar functions, preserves the 2d conformal gauge, since it 
affects the metric only by the change 



B(t,y)^B(i,y) 



B(t,y)-hn(f (t + y)g' (t - y)) 



(45) 

t,y-+t,y 



This is most easily seen in null coordinates, where 



v + u v — u 2 

t = — -=- , y = — -=- => as = —2 du dv . (46) 
V 2 y2 

Generally, the two branes will have a time dependent position in the extra dimension, 
described by the two functions y\ (t) and y 2 (t) , respectively. However, as long as their 
motion occurs slower than the speed of light, \yi\ and \y 2 \ < 1 , we can perform a change of 



coordinates (}44j) to have them at fixed position along y , as we now show. 

Let us first fix the first brane at y = . For this to happen, the two functions / and g 
appearing in equations (}4"4"|) have to satisfy 

f(t + y 1 (t))=g(t-y 1 (t)) . (47) 

We can choose / arbitrarily, and use (}4Tj) to determine g . The condition \yi\ < 1 guarantees 
this can be always done, since the arguments of both the two functions increase monotonically 
in time. 

In the new coordinate system, the first brane is fixed at y± = , while again the second 
one will be generally moving according to some function y 2 (t) . This function describes the 
parallel motion of the second brane with respect to the first one. Since in the old system 
of coordinates the relative motion was at a speed lower than that of light, this will be the 
case also in the new coordinates, | ^2 1 < 1 . To preserve the first brane at the origin, the 
residual freedom is restricted to / (w) = g (w) , i.e. / and g are the same functions of 
their arguments. If we choose / to satisfy 

f(t + y 2 (t)) = f(t-y 2 (t)) + 2, (48) 

we finally reach a third system of coordinates where the two branes are fixed at y% = and 
2/2 = 1, respectively. As before, the function / can always be constructed. Since \y 2 \ < 1 , 
the arguments of both terms increase monotonically in time. We can then use the value of 
/ at the right hand side of equation to "construct" the value of / at the left hand side. 



A. 2 Residual Gauge Freedom 

Even with the position of the branes fixed, the freedom of reparametrization is not exhausted 
yet. The residual gauge degrees of freedom are again of the form ([44)1 . with 

f(w) = g (w) =F(w) , F (w + 2) - F (w) = 2 . (49) 
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The most generic function F with this property is 

oo 

F(w) — w + (a n cos nxw + fc n sm nirw) , (50) 

n=0 

with arbitrary coefficients a n ,b n . The appearance of these gauge degrees of freedom is 
manifest in some of the numerical results we obtained, for example in Figure |U they are 
shown as ripples in the metric component B(t,y). 



A. 3 Perturbations of the Randall- Sundrum geometry 



It is interesting to note that, apart from pure gauge modes described in the previous ap- 
pendix, there exists only one kind of x— independent small fluctuations about the Randall- 
Sundrum geometry (|22j) . As we show now, this perturbation is related to a small change of 
the interbrane distance D , which is not stabilized without a scalar field. We know that any 
x— independent configuration can be written in the conformal gauge with the position of the 
branes at y = 0, 1 . Thus, all the perturbations we are interested in can be written in terms 
of the metric components A(t,y) = A (y) + 5a(t,y) , B(t,y) = B (y) + 5b(t,y) , where 
A = B is the Randall-Sundrum solution (|22|). To find which perturbations are allowed, 
we linearize the bulk Einstein equations. The dynamical ones reduce to 



8a — 5a" H h 



5b - 5b" 



y + i 

65a' 



(y + iY 

45b 



y + 7 (y + iY 



■■ o 





while the two constraint equations are conveniently recast in the form 



d ( , 5b 

— I 8a' + 

at V Z/ + 7 



d_ 

dy 



{y + 7) 3 5a' + 



5b 



y + i 



o , 

. 



Here 7 = (e D ^ — l) 1 , cf. equation (}2*2*j) . The last two equations give 



5a' + 



8b 



= C(y + 1 f 



y + i 

with C constant. Finally, linearizing the boundary conditions we have 

5b 



5a' + 



y + i 



, 



\y=o,i 



+ 



y + 7 







\y=o,i 



(51) 



(52) 



(53) 



(54) 
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the first of which enforces C — . Substituting equation (}53|) into the second equation of (}5T| . 
we have a differential equation in terms of 5b and its derivatives only. Fourier transforming 



5b 



due luJt 5b (u,y) 



we get an ordinary differential equation which is solved by 



5b = F 



COS (uj z) 



sm l u z) 



UJ z 



+ G 



sin {uj z) + 



cos [uj z) 



uj z 



z = y + >y. 



(55) 



(56) 



The boundary conditions for 5b at the two branes become two equations for the two 
parameters F and G . Non-vanishing solutions are possible only for 



uj = n it 



n = 0,1,2, 



(57) 



For these values, the two coefficients F and G are related by G = Ftanwy . The Fourier 
transform of 5a is then easily obtained from the remaining two equations 



5a = — [-F sin {uj z) + G cos {ujz)] + K5 {uj) 

UJ z 

where K is a constant. Back in coordinate space 

— F n e innt s in (n 7t y) 



(58) 



5a 
5b 



' cos {mvy) nrr (y + 7) 



+ K 



Annt 



n cos (rnij) 



cos (nny) — 



sin (niry) 
nn (y + 7) 



(59) 



These are the most generic x— independent perturbations of the Randall-Sundrum so- 
lution (|2*2*jl . However, most of them are pure gauge modes. Let us consider infinitesimal 
change of coordinates of the residual gauge discussed in Appendix IA.2I 

F(w) = w + F(w) = w + J2fne innw , 

n 

t -> t + J2fne innt cos(n7ry) , 

n 

V 2/ + ^/ne l ^'zsin(n7ry) . (60) 

n 

Under this infinitesimal change of coordinates, the metric coefficients undergo the infinites- 
imal changes 



A 
B 



F(t + y)-F(t-y) 

n T A-q , 

_ F(t + y)-F(t-y) P (t + y) + P (t - y) 

o o 



(61) 
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By choosing 

U = r , n = l,2,..., (62) 

n ir cos (nir'y) 

we see that the perturbations (JoTJj) are equivalent to 

A + 5a = A -^- + K = A + ^^ + K 

y + i y + i 

B Q + 5b = B + ^-, (63) 

y + i 

where /C = K — Fq is also constant. 

By an appropriate rescaling of the spatial x coordinates we can set K to zero. 



B Determination of Static Configurations 

Here we describe the method used to determine static solutions, once the bulk and brane 
potentials for the scalar field are given. This is done by a numerical boundary value problem 
solver using the shooting method. As discussed in Section 13.31 we first deal with the two 
boundary conditions at the first brane. Any two of B , <ft , <ft' Q , and H can be chosen freely, 
while the other two are determined by the junction conditions at the first brane. We find 
it more convenient to choose the values of B Q and 0' o , since the latter cannot be taken 
arbitrarily large if we want the solutions to remain regular all across the bulk. A fourth 
order Runge-Kutta integrator is then employed to integrate equations (jl9)l in the bulk. The 
aim is to find the values of Bq and <ft' for which the solutions are regular in the < y < 1 
interval, and for which the boundary conditions on the second brane are also satisfied. We 
can recast the latter in the form 



Cl = B[-±U 1 e* = , c 2 ^' 1 + ^e Bl =0. (64) 
6 2 d(p 

Both c\ and c 2 are only functions of the chosen value for Bq and <p' Q , and in general do 
not vanish. We use Newton's method to find the zeros of these two functions, that is the 
initial configurations at the first brane for which the junction conditions at the second brane 
are also satisfied. In practice, for the potentials we have studied, Newton's method does not 
converge globally. Fortunately, the bulk equations (fTT?j) can be integrated very quickly, so 
that we can perform a "brute force" scan in the {B , <p' } plane. We then apply Newton's 
method starting only from those values which are sufficiently close to a solution, i.e. for 
which ci and c 2 turn out to be sufficiently close to zero. 

The existence of static solutions is not guaranteed for arbitrary bulk and brane potentials. 
As discussed in [Ej, many potentials do not give static solutions at all, while some others 
typically lead to a finite number of them. Using the geometrical method of phase portraits 
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Figure 10: Phasespace illustration of the two solutions for given potentials 

for quadratic potentials (j2J) we found at most two static solutions in the (wide) space of 
possible initial configurations we have scanned. Figure shows how the phase portrait 
method allows us to visualize the quest for static configurations. Following the method 
of [T7j, we draw curves in the {<f>, e~ B (f>', e~ B B'} phase space. 

Each of the two thick curves refers to one of the two branes, and it joins points for which 
the junction conditions on that brane are satisfied. The thin curves are are a sampling of bulk 
trajectories which satisfy the junction conditions at one of the two branes (at y — 0). Valid 
static solutions consist of trajectories that satisfy the junction conditions at both branes. 
Hence, in the phase portrait they are represented by the few trajectories which intersect 
both the thick curves. Lines perpendicular to the trajectories represent lines of y = const. 
We see that in the case at hand, corresponding to quadratic potentials ©, there are two 
intersections. 

B.l Perturbations of Static Solutions 

Generic perturbations around a static configuration (fTSj) are described by the functions 
5B (y), SA (y), S(j) (y) and their first time derivatives 5B (y), 5A (y), and <50 o (?/); which 
are obtained by equating the time dependent fields and the first time derivatives at an 
initial moment to = 0. In the bulk, four of them can be specified arbitrarily, while the 
remaining functions are obtained from the constraint equations. One possible choice of 
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initial perturbations, adopted in the example of Figure is given by 

5B (y) = 6A o (y)=6A o (y)=0 , 

8<j) (y) = 5(f) (y) , 

Wo(y) = §0'- (^<f)' 2 + 5<f)' 2 -2<f)'5<f)'-2e 2B [V(<f))-V(<f)-5<f))] 

SB (y) = (65) 

The perturbation 5<f> (y) can be specified arbitrarily. In the example shown, the Gaussian 
profile 



(y - 1/2) 



2 



i —) , (66) 

is centered between the branes. A sufficiently small values for h guarantee that the pertur- 
bations are exponentially suppressed at the brane locations, leaving the junction conditions 
(practically) unaffected. 
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Table 1: Parameters used for simulations presented in the figures 
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